/* PROGRAM EDGE

   --- edge [compcrit] < input > output
   --- computes the distance to the nearest edge for each galaxy
      *[compcrit] : critical completeness that defines survey boundary 
                   (default=0.9)
      * < lss sample (e.g., lss.sample10bright13.dat)
      * > output : ascii file containing <ID RA DEC cz rangle rdist>
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "mask.h"
#define pi 3.1415926535898
#define sqr(x) ((x)*(x))

main(int argc, char *argv[])
{
  int i,j,junki,ID ;
  double junkf,RA,DEC,cz ;
  double theta,phi,x[3],comp,compcrit ;
  double sintheta,costheta,sinphi,cosphi ;
  double r,dr,rmin,rmax,sinr,cosr,psi,dpsi,sinpsi,cospsi ;
  double sintheta2,costheta2,theta2,cosDphi,Dphi,phi2 ;
  double comp2,redge ;
  char *maskfile;
  /*double completeness(double *) ;*/

/*---Read-Arguments---------------------------------------------------*/
  maskfile = argv[1];

  compcrit=0.9 ;
  if(argc>2) sscanf(argv[2],"%f",&compcrit) ;

/*---Load-Mask-Info---------------------------------------------------*/
  loadmask(maskfile) ;

/*---Set-Grid-Parameters----------------------------------------------*/
  rmin = 1.0 ;        /*---In-arcminutes---*/
  rmax = 200.0 ;
  dr = 1.0 ;
  dpsi = 5.0 ;        /*---In-degrees---*/

  rmin *= pi/(60*180.0) ;  /*---Convert-to-radians---*/
  rmax *= pi/(60*180.0) ;
  dr *= pi/(60*180.0) ;
  dpsi *= pi/180.0 ;

/*---Read-Galaxy-file-------------------------------------------------*/
  i = 0 ;
  while(fscanf(stdin,"%d %d %d %f %f %f %f %f",&ID,&junki,&junki,&RA,&DEC,&cz,&junkf,&junkf)!=EOF)
    {
      fprintf(stderr,"Doing %d\n",ID) ;
      theta = (90-DEC)*pi/180.0 ;
      phi = RA*pi/180.0 ;

      x[0] = sin(theta)*cos(phi) ;
      x[1] = sin(theta)*sin(phi) ;
      x[2] = cos(theta) ;
 
      comp = completeness(x) ;

      if(comp>compcrit)
	{
	  sintheta = sin(theta) ;
	  costheta = cos(theta) ;
	  sinphi = sin(phi) ;
	  cosphi = cos(phi) ;

	  comp2 = 1.0 ;
	  for(r=rmin;r<rmax;r+=dr)
	    {
	      sinr = sin(r) ;
	      cosr = cos(r) ;

	      for(psi=0;psi<pi;psi+=dpsi)
		{
		  cospsi = cos(psi) ;

		  costheta2 = costheta*cosr + sintheta*sinr*cospsi ;
		  theta2 = acos(costheta2) ;
		  sintheta2 = sin(theta2) ;

		  cosDphi = (cosr - costheta*costheta2)/(sintheta*sintheta2) ;
		  Dphi = acos(cosDphi) ;
		  phi2 = phi + Dphi ;

		  x[0] = sintheta2*cos(phi2) ;
		  x[1] = sintheta2*sin(phi2) ;
		  x[2] = costheta2 ;

		  comp2 = completeness(x) ;

		  if(comp2<compcrit) 
		    break ;

/*---Check-opposite-side-of-circle:-180-psi---*/		  
		  phi2 = phi - Dphi ;
		  x[0] = sintheta2*cos(phi2) ;
		  x[1] = sintheta2*sin(phi2) ;
		  x[2] = costheta2 ;
		  
		  comp2 = completeness(x) ;

		  if(comp2<compcrit) 
		    break ;
		}
	      if(comp2<compcrit) 
		break ;
	    }
	  r = 180*60*r/pi ;
	  redge = cz*sinr ;
	}
      else
	{
	  r = -9.0 ;
	  redge = -9.0 ;
	}
      fprintf(stdout,"%7d %10.6f %10.6f %8.2f %8.6f  %8.2f %8.2f\n",ID,RA,DEC,cz,comp,r,redge) ;
    }
}










